Strong-coupling effects in the relaxation dynamics 
of ultracold neutral plasmas 



T. Pohl and T. Pattard 
February 2, 2008 



Max Planck Institute for the Physics of Complex Systems, 
Ndthnitzer Str. 38, D-01187 Dresden, Germany 

Abstract 

We describe a hybrid molecular dynamics approach for the description of 
ultracold neutral plasmas, based on an adiabatic treatment of the electron 
gas and a full molecular dynamics simulation of the ions, which allows us to 
follow the long-time evolution of the plasma including the effect of the strongly 
coupled ion motion. The plasma shows a rather complex relaxation behavior, 
connected with temporal as well as spatial oscillations of the ion temperature. 
Furthermore, additional laser cooling of the ions during the plasma evolution 
drastically modifies the expansion dynamics, so that crystallization of the ion 
component can occur in this nonequilibrium system, leading to lattice-like 
structures or even long-range order resulting in concentric shells. 

1 Introduction 

Experiments in cooling and trapping of neutral gases have paved the way toward 
a new parameter regime of ionized gases, namely the regime of ultracold neutral 
plasmas (UNPs). Experimentally, UNPs are produced by photoionizing a cloud of 
laser-cooled atoms collected in a magneto-optical trap PP , with temperatures down 
to 10 /iK. By tuning the frequency of the ionizing laser, initial electron kinetic 
energies of E c /k-Q = IK — 1000K have been achieved. The time evolution of several 
quantities characterizing the state of the plasma, such as the plasma density ^EUHHi 
the rate of expansion of the plasma cloud into the surrounding vacuum [3], the 
energy-resolved population of bound Rydberg states formed through recombination 
0], or electronic |S| as well as ionic |3] temperature have been measured using 
various plasma diagnostic methods. 

Despite the low typical densities of » 10 9 cm -3 , the very low initial temperatures 
suggest that these plasmas have been produced well within the strongly coupled 
regime, with Coulomb coupling parameters up to T c = 10 for the electrons and 
even T; = 30000 for the ions. Thus, UNPs seem to offer a unique opportunity for 
a laboratory study of neutral plasmas where, depending on the initial electronic 
kinetic energy, either one component (namely the ions) or both components (ions 
and electrons) may be strongly coupled. Moreover, the plasma is created in a 
completely uncorrelated state, i.e. far away from thermodynamical equilibrium. The 
relaxation of a strongly correlated system towards equilibrium is an interesting topic 
in non-equilibrium thermodynamics and has been studied for decades. The history 
of this problem must be traced back to the important contributions of Klimontovich 
HI 00 IHi wn0 pointed out that kinetic energy conserving collision integrals such as 
the Boltzmann, Landau and Lenard-Balescu collision integrals are not appropriate 
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for such a situation, and derived non-Markovian kinetic equations taking correctly 
into account total energy conservation of the system. In the following years this 
problem has attracted much attention and the relaxation of nonequilibrium strongly 
coupled plasmas has been studied by a variety of different methods ^3 ^2 
I13| . The very low densities of UNPs make it now possible to directly observe 
the dynamical development of spatial correlations, which may serve as the first 
experimental check of the present understanding of the strongly coupled plasma 
dynamics. Moreover, it turns out that the timescale of the plasma expansion, 
the correlation time as well as the relaxation time of the ions are almost equal. 
Therefore Bogoliubov's functional hypothesis, usually used in kinetic theory, breaks 
down under the present conditions, which may lead to a very interesting relaxation 
behavior but also causes some difficulties in the theoretical description of these 
systems, since the plasma dynamics can not be divided into different relaxation 
stages. 

2 Theoretical approach 

A full molecular dynamics simulation of ultracold neutral plasmas over experimen- 
tally relevant timescales is infeasible with present-day computer resources due to 
the large number of particles (N « 10 5 ) and the long observation times (t w 10~ 4 
s) involved. In order to model the evolution of UNPs, we have developed a hybrid 
molecular dynamics (HMD) approach which treats electrons and ions on different 
levels of sophistication, namely in a hydrodynamical approximation on the one hand 
(for the electrons) and on a full molecular dynamics level on the other hand (for the 
ions) [T3|. For the electrons, it has been shown that several heating effects, such 
as continuum threshold lowering |15|. build-up of correlations |16|. and, predomi- 
nantly, three-body recombination |17| rapidly increase the electronic temperature. 
As a consequence, the electrons are always weakly coupled, r c < 0.2, over the 
whole course of the system evolution. Moreover, due to the small electron-to-ion 
mass ratio, the relaxation timescale of the electrons is much smaller than that of 
the ions as well as the timescale of the plasma expansion. Hence, an adiabatic ap- 
proximation may safely be applied, assuming instant equilibration of the electrons. 
This allows for the use of much larger timesteps than in a full MD simulation since 
the electronic motion does not need to be resolved. It is this adiabatic approxi- 
mation for the electrons which makes a molecular dynamics treatment of the ionic 
motion in UNPs computationally feasible. The electronic density is determined self- 
consistently from the Poisson equation. The fact that the potential well created by 
the ions which is trapping the electrons has a finite depth is taken into account 
by using a King-type distribution |18j known from simulations of globular clusters 
rather than a Maxwell-Boltzmann distribution for the electron velocities, with the 
electronic temperature T e obtained from energy conservation. The finite well depth 
also leads to evaporation of a fraction of the free electrons in the very early stage of 
the plasma evolution, which is accounted for by determining the fraction of trapped 
electrons from the results of pQ . The dynamics of the heavy particles is described 
in the framework of a chemical picture, where inelastic processes, namely three- 
body recombination and electron impact ionization, excitation and deexcitation, 
are taken into account on the basis of Boltzmann-type collision integrals |191 |S] , 
with the transition rates taken from [201 ■ Numerically, the resulting collision in- 
tegrals are evaluated using a Monte Carlo sampling as described in El 122] • 
The ions and recombined atoms are then propagated individually in a molecular 
dynamics simulation, taking into account the electronic mean-field potential and 



the full interaction potential of the remaining ions 1 . In order to allow for larger 
particle numbers, the most time-consuming part of the HMD simulation, namely 
the calculation of the interionic forces, is done using a treecode procedure originally 
designed for astrophysical problems , which scales like Ni In Ni rather than Nf 
with the number N of ions. 

As shown in several publications |14l 1241 TIFH I22| . the HMD approach outlined 
above provides a powerful method for the description of UNPs, taking full account 
of ionic correlation effects. However, due to the large numerical effort involved, it 
is limited to particle numbers of N ~ 10 5 . While this permits a direct simulation 
of many, particularly of the early, experiments, an increasing number of experi- 
ments is performed with larger particle numbers up to 10 7 . Thus, an alternative 
method which is able to treat such larger systems is desirable. Such a method is 
indeed available ^1], based on a hydrodynamical description of both electrons and 
ions similar to that introduced in ^| |2] ■ Starting from the first equation of the 
BBGKY hierarchy, one obtains the evolution equations for the one-particle distri- 
bution functions / of the electrons and ions. Neglecting again electron-electron as 
well as electron-ion correlations, and employing the same adiabatic approximation 
for the electrons already used in the HMD approach, a quasineutral approximation 
[2l)| permits expressing the mean-field electrostatic potential in terms of the ionic 
density, leading to a closed equation for the ion distribution function which contains 
the electron temperature as a parameter. A Gaussian ansatz for the ion distribution 
function, 



which corresponds to the initial state of the plasma cloud, is then inserted into 
the evolution equations for the second moments (r 2 ), (rv) and (v 2 ) of the ion dis- 
tribution function. In this way, evolution equations for the width a of the cloud, 
the parameter 7 of the hydrodynamical expansion velocity u = ^yr and the ionic 
temperature 7] are obtained. Ionic correlations are taken into account in an approx- 
imate way using a local density approximation together with a gradient expansion, 
reducing the description of their influence on the plasma dynamics to the evolution 
of a single macroscopic quantity, namely the correlation energy Ua of a homoge- 
neous plasma. The relaxation behavior of Ua is modeled using a correlation-time 
approximation 27 with a correlation time equal to the inverse of the ionic plasma 
frequency, T cor r = j\ together with an analytical expression for the equilibrium 
value of Ua [2E1- Finally, inelastic processes such as three-body recombination and 
electron impact ionization, excitation and deexcitation are incorporated on the ba- 
sis of rate equations, and the influence of the recombined Rydberg atoms on the 
expansion dynamics is taken into account assuming equal hydrodynamical velocities 

1 In order to bring out clearly the role of ionic correlations, it is also possible to neglect them 
in the HMD approach by propagating the ions in the mean-field potential created by all charges 
rather than the full ionic interaction. 
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Figure 1: Electronic temperature T e (t) for an expanding plasma of 40000 Sr ions 
with an initial average density of pi — 10 9 cm -3 and an initial electron kinetic energy 
of 20 K, obtained from the HMD simulation (a) and from eqs. (j2J) (b), with (solid) 
and without (dotted) the inclusion of ionic correlations. The inset shows the ratio 
of the electron temperatures obtained from both methods. 



for atoms and ions. The final set of evolution equations then reads 

a = 7 cr , (2a) 
Ni (k B T e + k B Ti + ±U a ) 2 

7 = 7 ' (2b) 

k B t = -2jk B Ti - pUa - , (2c) 

Ua = -w P;i (c7H-^ eq) ) (2d) 
■Mi(n) = [R hh (p,n)Af a (p) ~ R hh {n,p)Af a (n)} 

+i?tbr(n)iV i - i? ion (n)A/;(n) (2e) 
and the electronic temperature is determined by energy conservation, 

N,k B T c + [Ni + N a ] [k B T { + m n 2 a 2 ] + -Ntf* - H y*JV a (n)^ = const. , (2f) 

n 

where Af a (n) defines the population of Rydberg states, iV a = JZ n is the total 

number of atoms and 1Z = 13.6eV is the Rydberg constant. 

The preceeding hydrodynamical method is much more approximate than the 
HMD approach, but, on the other hand, it is much simpler and quicker. For parti- 
cle numbers of N t « 10 5 , it requires about two orders of magnitude less CPU time. 
Since its computational effort is independent of the number of particles, it allows for 
a simulation of larger plasma clouds corresponding to a number of current exper- 
iments. Moreover, and maybe equally important, it provides physical insight into 
the plasma dynamics since it is based on a few simple evolution equations for the 
macroscopic observables characterizing the state of the plasma. As we have inves- 
tigated in detail in ^3] > there is generally surprisingly good agreement between the 
hydrodynamical simulation and the more sophisticated HMD calculation as long 
as macroscopic, i.e. spatially averaged, quantities such as electronic temperature, 
expansion velocity, ionic correlation energy etc. are considered. As an example, we 
show in figure ^ the time evolution of the electronic temperature for a plasma of 
40000 Sr ions with an initial average density of 10 9 cm~ 3 and an initial electron 
kinetic energy of 20 K, obtained from the HMD simulation (a) and from eqs. (J2J 
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Figure 2: Electronic temperature T e (t) for a plasma of 1.2 • 10 6 Xenon ions with 
an initial average density of 1.35 • 10 9 cm -3 for two different initial electron temper- 
atures, T c = 6.67K (filled dots) and T c = 66.67K (open dots). The lines show the 
hydrodynamical simulation (solid lines: including inelastic collisions, dashed lines: 
without inelastic collisions), the dots the experiment [H], scaled down by 26% (see 
text). 

(b). During the whole system evolution, the agreement between the two simulation 
methods is better than about 8%, and it becomes even better at later times. Thus, 
we conclude that, for the present type of experimental setups, the hydrodynami- 
cal method outlined above, and in particular the approximate treatment of ionic 
correlations, is well suited for the description of the behavior of UNPs. 

3 Results and discussion 
3.1 Comparison with experiments 

In fact, fig. n only shows good agreement between the two theoretical simulation 
methods, without comparison with experiment. Such a comparison is now also 
possible, since measurements of the electron temperature dynamics have recently 
been reported in [5]- Fig. [2] shows the time evolution of the electronic tempera- 
ture for a Xenon plasma with Ni(0) = 1.2 • 10 6 , pi(0) — 1.35 • 10 9 cm~ 3 and two 
different initial temperatures of T e (0) = 66.67K and T c (0) = 6.67K. In addition 
to the full hydrodynamical simulation according to equations J5J, fig. also shows 
corresponding calculations where the effect of inelastic electron-ion collisions, eq. 
(j2"e)l . is neglected (dashed lines). (The plasmas in these experiments are too large 
to be simulated using the HMD approach.) For the high initial temperature, 
there is close agreement between the two corresponding simulations, showing that 
inelastic processes are almost negligible in this case. Indeed, it is known that the 
high-temperature plasma expansion is well described by the collisionless plasma 
dynamics, and the hydrodynamical model is expected to accurately reproduce the 
plasma dynamics in this regime. Since an overall systematic error of about 70% for 
the temperature measurement has been reported in 0, we have exploited this fact 
to calibrate the measured temperatures by scaling down both experimental data 
sets by 26% in order to match the high-temperature results to our calculations. As 
can be seen in the figure, there is excellent agreement between simulation and ex- 
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Figure 3: Time evolution of the average electron density of a Xenon plasma of 
500000 ions with an initial average density of 10 9 cm~ 3 and an initial electron tem- 
perature of T c = 210K (a) and T c — 2.6K (b). The lines show the results of the 
model equations (solid lines: including inelastic collisions, dashed lines: without 
inelastic collisions) and the dots the experimental data from 0. 

periment also for the lower temperature. (We stress that there is no further scaling 
of the low-temperature data in order to achieve quantitative agreement, the same 
calibration factor as in the high-temperature case is used.) In this case, inelastic 
collisions play a decisive role for the evolution of the system. More specifically, 
as has been found already in |17j . three-body recombination heats the plasma and 
significantly changes its behavior, leading to a weakly coupled electron gas, as dis- 
cussed above in connection with the omission of electronic correlation effects in the 
numerical treatment. Moreover, there has been some discussion in the literature 
whether the collision rates of [30] would still be applicable at these ultralow tem- 
peratures, or whether three-body recombination would be significantly altered. The 
close agreement between the present simulation and the experimental data in fig. El 
suggests that the rates of [20], while ultimately diverging oc T<r 9 ^ 2 for T c — > 0, still 
adequately describe three-body recombination processes in the temperature range 
under consideration. 

As a second example, figure E3 shows the time evolution of the electronic density 
for a Xenon plasma of 500000 ions with an initial average density of 10 9 cm~ 3 and two 
different initial electron temperatures of T o (0) = 210K and T e (0) = 2.6K gj. Again, 
it can be seen that the model equations nicely reproduce the density evolution in 
both temperature regimes, in agreement with ^7] where it was shown that the low- 
temperature enhancement of the expansion velocity 2 is caused by recombination 
heating and is not due to strong-coupling effects of the electrons. 

3.2 Role of ionic correlations 

Having thus established the validity of our numerical methods for the description 
of UNPs, we can now turn to a more detailed investigation of the role of ionic cor- 
relations in these systems. It is found that, for situations corresponding to the type 
of experiments EJ Ej , they hardly influence the macroscopic expansion behavior 
of the plasma. This becomes evident, e.g., in fig.H where the "full" simulations as 
described above (solid lines) are compared to a mean-field treatment of the system 
completely neglecting correlation effects (dotted lines). The correlation- induced 
heating of the ions |30] 1301 1ST] leads to a slightly faster expansion of the plasma, 
which in turn results in a slightly faster adiabatic cooling of the electrons |14|. 
However, the overall effect is almost negligible. 

A closer look, on the other hand, reveals that certain aspects of the expansion 
dynamics are indeed significantly affected by the strong ion-ion interaction, as can 
be seen in figure 0] There, the spatial density of the ions is shown after t = 3 /zs 
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Figure 4: Spatial density p\ (solid) of the ions, at t = 3 fis, compared to the 
Gaussian profile assumed for the kinetic model (dashed). Additionally, p\ obtained 
from the particle simulation using the mean-field interaction only is shown as the 
dotted line. Initial-state parameters are the same as in fig. ^ 

for the same plasma as in fig. ^ A mean-field treatment of the particle interactions 
PTj predicts that a shock front should form at the plasma edge, seen as the sharp 
spike in fig. 0] (dotted line). Apparently, with ionic correlations included (solid line) 
the peak structure is much less pronounced than in mean-field approximation. This 
is due to dissipation caused by ion-ion collisions which are fully taken into account 
in the HMD simulation. As shown in |.'S2|. by adding an ion viscosity term to the 
hydrodynamic equations of motion, dissipation tends to stabilize the ion density 
and prevents the occurrence of wavebreaking which was found to be responsible for 
the diverging ion density at the plasma edge in the case of a dissipationless plasma 
expansion. Furthermore, the initial correlation heating of the ions largely increases 
the thermal ion velocities, leading to a broadening of the peak structure compared 
to the zero-temperature case. 

Another obvious aspect where ionic correlations play a dominant role is the 
behavior of the ionic temperature. Considering the huge ionic coupling constants 
suggested by the low initial ion temperatures, this temperature turns out to be 
an important quantity since it directly determines the value of Y\. According to a 
mean-field treatment, the ions would remain the (near) zero temperature fluid they 
are initially. However, as has been pointed out before, the ions are created in a com- 
pletely uncorrelated non-equilibrium state, and they quickly heat up through the 
build-up of correlations as the system relaxes toward thermodynamical equilibrium. 
As shown in |14| . even at early times the ionic velocity distribution is locally well 
described by a Maxwell distribution corresponding to a (spatially) local tempera- 
ture, justifying the definition of a — due to the spherical symmetry of the plasma 
- radius-dependent ion temperature Ti(r,t). Moreover, if the spatially averaged 
temperature is identified with the ion temperature determined by the model equa- 
tions one can find again good agreement between both approaches concerning 
the timescale of the initial heating as well as the magnitude of the ion temperature, 
even at later times However, as becomes apparent from fig. [3 the HMD simula- 
tions show temporal oscillations of the ionic temperature, which can, of course, not 
be described by the linear ansatz of the correlation-time approximation. Such tem- 
poral oscillations of the temperature during the initial relaxation stage are known 
from molecular dynamics simulations of homogeneous one-component and two- 
component plasmas, which are clearly caused by the strongly coupled collective 




Figure 5: Time evolution of the density-scaled average ionic temperature for a 
plasma consisting of 400000 ions with an initial electronic Coulomb coupling pa- 
rameter of T c (0) = 0.07. 




Figure 6: Time evolution of the ion temperature determined from a central sphere 
with a radius of half of the plasma width a (a) and time dependence of the amplitude 
of the corresponding oscillations (b). The initial-state parameters are the same as 
in fig. [3 

ion dynamics, since they increase in strength with increasing Tj and disappear for 
r;(0) < 0.5 Q31. 

Despite the fact that the maximum initial coupling constant used in |13| is 
Ti = 5, while a value of r;(0) « 40000 is considered in the case of fig. El the oscilla- 
tions observed in are much more pronounced and persist much longer than in 
the present case. It becomes apparent that the rapid damping of the temperature 
oscillations can be traced to the inhomogeneity of the Gaussian density profile by 
looking at the central part of the plasma only, where the ionic density is approxi- 
mately constant (figure |BJ). The temperature oscillations with an oscillation period 
of half of the inverse plasma frequency vo = uj p ^(0)/2tt defined in the central plasma 
region are much more pronounced in this case, showing an exponential decay with a 
characteristic damping rate of vq (fig. E^b)). The temporally and spatially resolved 
temperature evolution shown in fig.0shows that the radially decreasing ion density 
leads to local temperature oscillations with radially increasing frequencies, thereby 
causing also spatial oscillations of the local ion temperature. Therefore, the seem- 
ingly enhanced damping rate, which has also been observed in recent experiments, 
is purely an effect of the averaging of these local oscillations over the total plasma 



Figure 7: Temporally and spatially resolved time evolution of the ion temperature. 
The initial-state parameters are the same as in fig. [5] 

volume. 

3.3 Coulomb crystallization through laser cooling 

The above considerations show that, while not dramatically affecting the overall 
expansion behavior of the plasma cloud, strong-coupling effects play an important 
role in different aspects of the evolution of UNPs. Thus, UNPs provide a prime ex- 
ample of laboratory realizations of strongly nonideal plasmas. Moreover, the HMD 
approach developed in |14| is well suited for an accurate description of these systems 
over experimentally relevant timescales, allowing for direct comparison between ex- 
periment and theory. Many interesting aspects of the relaxation behavior of these 
non-equilibrium plasmas may thus be studied in great detail. However, while effects 
of strong ionic coupling become apparent in UNPs, the naively expected regime with 
r > 100 can not be reached with the current experimental setups. For the elec- 
trons, it is predominantly three-body recombination which heats them by several 
orders of magnitude, so that T c < 0.2 during the whole system evolution. The ionic 
component, on the other hand, is heated by the correlation- induced heating until 
T; w 1, i.e. just at the border of the strongly coupled regime I2HH2]- Thus, it is the 
very build-up of ionic correlations one wishes to study that eventually shuts off the 
process and limits the amount of coupling achievable in these systems. 

As soon as the reason for this ionic heating became clear, several proposals have 
been made in order to avoid or at least reduce the effect, among them (i) using 
fcrmionic atoms cooled below the Fermi temperature in the initial state, so that the 
Fermi hole around each atom prevents the occurrence of small interatomic distances 
|29| : (ii) an intermediate step of exciting atoms into high Rydberg states, so that 
the interatomic spacing is at least twice the radius of the corresponding Rydberg 
state |80| : and (Hi) the continuous laser-cooling of the plasma ions after their initial 
creation, so that the correlation heating is counterbalanced by the external cooling 
(341 124| . We have simulated the latter scenario using the HMD method, extended 
to allow for the description of laser cooling, as well as elastic electron-ion collisions 
which are negligible for the free plasma expansion but not necessarily in the laser- 
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Figure 8: Radial density and a central slice of a plasma with Ni(0) — 80000, 
T e (0) = 0.08, cooled with a damping rate of f3 = 0.2oj Pi i(0) at a time of w Pj ;(0)i = 
216. (For better contrast, different cuts have been overlayed.) 




Figure 9: Arrangement of the ions on the first (a), third (b) and fifth (c) shell of 
the plasma of fig. |H1 

cooled case [21 122] ■ Laser cooling is modeled by adding a Langevin force 

Fcooi = -rm0v + v / 2pk B T c m4 (3) 

to the ion equation of motion, where v is the ion velocity, £ is a stochastic variable 
with (£) = 0, (£(t)£(t + T)) = 3<5(r), and the cooling rate {3 and the correspond- 
ing Doppler temperature T c are determined by the properties of the cooling laser 
[3"K] . Elastic electron-ion collisions are taken into account on the basis of the corre- 
sponding Boltzmann collision integral, which is again evaluated by a Monte-Carlo 
procedure [22J - 

It is found that laser cooling leads to qualitative changes of the plasma dynam- 
ics. In particular, it significantly decelerates the expansion of the plasma, whose 
width is found to increase only as a oc i 1 ' 4 , in contrast to freely expanding plasmas 
which behave as a oc t. It is this drastic slow-down of the expansion which favors 
the development of strong ion correlations, compared to a free plasma where the 
expansion considerably disturbs the relaxation of the system. The simulations show 
further that strongly coupled expanding plasmas can indeed be created under real- 
istic conditions, with ionic coupling constants far above the crystallization limit for 
homogeneous plasmas of Ti « 174 PS). Here we find, depending on the initial con- 
ditions, i.e. ion number and initial electronic Coulomb coupling parameter, strong 
liquid-like short-range correlations or even the onset of a radial crystallization of 
the ions. This is demonstrated in fig. |S| showing the radial density and a central 
slice of a plasma with Ni(0) = 80000, r e (0) = 0.08, cooled with a damping rate 
of P = 0.2w P! i(0), at a scaled time of u) p ^(0)t — 216. The formation of concentric 



shells in the center of the cloud is clearly visible. As illustrated in fig. El beside the 
radial ordering there is also significant intra-shell ordering, namely a formation of 
hexagonal structures on the shells, which are, however, considerably disturbed by 
the curvature of the shells. 

4 Conclusions 

In summary, we have used an HMD approach to study the behavior of ultracold 
neutral plasmas on long time scales. We have shown that effects of strong interionic 
coupling are indeed visible in such systems, e.g. most prominently in the relaxation 
behavior of the ion temperature, which is connected with transient temporal as well 
as spatial oscillations. Nevertheless, the strongly coupled regime of T > 100 is not 
reached with the current experimental setups. We have demonstrated, however, 
that additional continuous laser cooling of the ions during the plasma evolution 
qualitatively changes the expansion behavior of the system and should allow for the 
Coulomb crystallization of the plasma |241 \22\ . It will be an interesting subject for 
further investigation to study in detail the dynamics of this crystallization process, 
which differs from the shell structure formation observed in trapped nonneutral plas- 
mas |36| as explained in |24) . In particular, the influence of the plasma expansion, 
which presumably causes the transition from liquid-like short-range correlation to 
the radial ordering, deserves more detailed studies. Other future directions include 
the study of effects induced by additional magnetic fields, or of ways to confine the 
plasma in a trap. 

We gratefully acknowledge many helpful discussions with J.M. Rost, as well as 
conversations with T.C. Killian and F. Robichcaux. 
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